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Abstract 

The article presents a series of numerical simulations of exact solutions of the 
Einstein equations performed using the Cactus code, a complete 3-dimensional ma- 
chinery for numerical relativity. We describe an application ("thorn") for the Cactus 
code that can be used for evolving a variety of exact solutions, with and without 
matter, including solutions used in modern cosmology for modeling the early stages 
of the universe. Our main purpose has been to test the Cactus code on these well- 
known examples, focusing mainly on the stability and convergence of the code. 



1 Introduction 



Numerical Relativity is concerned with the study of numerical solutions of the Einstein's 
equations for the gravitational field, which are at the core of the theory of General Rel- 
ativity. General Relativity is a 4-dimensional theory involving one dimension of 
time and three of space. The field equations, called Einstein equations, are : 
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where A is the cosmological constant, Rij the Ricci tensor, R the Ricci scalar, gij the 
spacetime metric, the stress-energy tensor, G the gravitational constant, c the speed 
of light and i,j = 0, 1, 2, 3. 

These equations are an extremely complicated system of coupled, non-linear, partial 
differential equations and solving them numerically makes enormous demands on the pro- 
cessing power and memory of a computer. Because of this. Numerical Relativity (0, 0) 
has proceeded in several stages, first solving 1-dimensional problems (that is 1 spatial di- 
mension, e.g. spherical symmetry), and then moving to 2-dimensional problems (i.e. axial 
symmetry). Only in the last few years have computers developed sufficiently to consider 
tackling fully 3-dimensional problems. At the Albert Einstein Institute, an international 
team led by Edward Seidel has developed a complete 3D code for numerical relativity, 
which has been named the "Cactus code" (see [^]). The Cactus code has been mainly 
designed as a computational toolkit (freely available for the scientific community) for 
simulating different systems of partial differential equations. In the particular of relativ- 
ity, the code can be used to simulate systems with strong gravitational fields: collapsing 
gravitational waves, colliding black-holes, neutron stars, and other violent astrophysical 
processes generating gravitational waves [0]-|]13|. 



Among the main problems related to the development of the Cactus code has been 
that of testing the code for stability and convergence during the simulations. Several 
applications of the code (named generically "thorns") were designed for this purpose, 
one of them being the so-called thorn "Exact", were some known exact solutions of the 
vacuum Einstein equations were implemented. The original version of this thorn was 
written and has been furthered developed by many people. Thorn Exact was designed for 
a comparative study of the numerical evolution of exact solutions of the vacuum Einstein 
equations 0. The thorn requires the full 4-dimensional metric of a given exact solution 
in the whole spacetime. Its routines then generate 3-1-1 data (i.e. lapse function, shift 
vector, spatial metric, and extrinsic curvature) from the given exact solution both as 
initial data for the numerical evolution, and also at every iteration step of the evolution 
so that it can be directly compared with the numerically evolved data. A limited number 
of exact solutions were included (among them Minkowski, Novikov and several Black 
Hole spacetimes), but since the thorn was constructed in very modular form it has been 
easy to add new exact solutions by just writing a special subroutine for each new case 
(with the corresponding parameters). Thus we have added a series of space-times which 
are of cosmological interest, such as De Sitter, Friedmann, Kasner, DeMilne and Godel 
models. The main problem was that some of these models are exact solutions of the 
Einstein equations with matter, so the original thorn has to be generalized to allow one 
to introduce the components of the corresponding stress-energy tensor. 
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The current version of this thorn Exact now contains also several space-time models 
used in cosmology for treating the early stages of the universe. Thus, through this thorn, 
the Cactus code can be used for numerical cosmology studies, an example being the 
inflationary cosmology where inflation is controlled by a scalar field (this is the main 
purpose of our future work). 

This article is a report on the current status of thorn Exact. We shall present some 
of the numerical results obtained on running the Cactus code for several cosmological 
models introduced in this thorn. We point out here also how the code produces conver- 
gent numerical results and how stable the evolution can be for different cases. At the 
moment the Cactus code can use two different formulations for numerical integrating the 
Einstein equations , namely the standard Arnowitt-Deser-Misner formulation (im- 
plemented as the ADM thorn) and a recent reformulation introduced by Baumgarte and 
Shapiro ([|ll|], H^) based on previous work of Shibata and Nakamura (1|15[) (implemented 
in the BSSN thorn). The comparative study of these two formulations shows the better 
stabihty properties of the BSSN formulation (as already pointed out in |11|). A special 
point to mention here is the problem of the boundary behavior of the numerical simu- 
lations. The boundary conditions currently available in the Cactus code are specifically 
adapted for the simulations of either asymptotically fiat spacetimes, or periodic space- 
times. Because of this we have been forced to use special boundary conditions in thorn 
Exact which force the code to take the boundary values from the exact solution itself, 
except in those cases were we can use the "fiat" boundary condition already implemented 
in the code (the fiat boundary condition means that the value of the given field or fields 
at the boundary is simply copied from the value one grid point in along the direction 
normal to the boundary). 

In the next sections we shall present several results obtained after running Cactus for 
different metrics given by thorn Exact, on different computer architectures. The simula- 
tions were performed mainly on single processor machines, using both the AEI computer 
network (SGI or Dec machines with UNIX operating system) and a Pentium III machine 
with a UNIX FreeBSD operating system at the West University of Timi§oara. Some of 
the simulations, with similar results, were also done on the Origin 2000 supercomputer 
at the AEI. As a conclusion of our work we have extended the Exact thorn to deal with 
exact solutions of the Einstein equations with matter, in particular several cosmological 
solutions. The thorn is available, via internet downloading, for anyone interested in using 
it. 

Through this article and in the Cactus code we shall use geometrical units with G = 
c= 1. 
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2 De Sitter cosmologies 



Modern cosmology is based on the Robertson-Walker metric ([^,|jT^) having the line 
element as : 



ds' = -dt^ + R{tf 



1 — kr"^ 



(2) 



in spherical coordinates {t,r,6,(f)). Here R{t), the scale factor, is depending only on 
the time t. The parameter k has three values, namely k = 0,1, —1 for fiat, closed or 
respectively open spatial three geometry. The Einstein equations (|l|) for A = and a 
stress-energy tensor as for a perfect fluid 

Tij = {p + p)uiUj +pgij 

(where p is the pressure, p the density and Ui the velocity in comoving coordinates) have, 
as an exact solution, for a matter dominated universe {p = 0) and fiat spatial geometry 
{k = 0) the DeSitter metric (also called Einstein-DeSitter metric, see |T|, 0, ||T6[ and 
[|T7|1 ) where R{t) ~ t^/^. Thus DeSitter metric has the line element as : 

ds^ = -df + at^/^ J^^2 ^ ^2 (^^q2 ^ sin(^)2^02^j _ 

where a is a constant. Here the only non- vanishing component of the the stress-energy 
tensor is : 

Next we shall present some of the results obtained after running Cactus for the De 
Sitter metric. First we show the behavior of the radial metric function [g^r = R{ty/{1 — 
kr"^) in the above line element, being denoted with grr at the Cactus output) obtained 
after 20 iterations on a grid of 28^ points, with < x,y, z < 1 and At = 0.25 Ax (Fig. |1|). 
Here and in the next figures, "normalized" means the L2 norm of the respective function 
calculated using all the values of the function on the computational grid at one time and 
the respective output Cactus files are denoted with _nm2 extension. 

Cactus code is evolving the Einstein equations using the 3+1 decomposition of space- 
time (both in ADM or BSSN evolution methods). Thus the Einstein equations can be split 
in two groups : the dynamic equations (for the time derivatives of the three-dimensional 
metric and extrinsic curvature) and the constraint equations (the Hamiltonian constraint 
and the Momentum constraint) - see and The constraint equations are satisfied 
during all time evolution of the system. Thus one of the main tests on the convergence in 
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DeSitter spacetime DeSltter spacetime 

Evolution of the grr values of the grr at different times on the grid 
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time Diagonal line of the grid 

Figure 1: Evolution of the L2 norm of the radial metric function g^r versus time (left 
panel) and of the metric function itself through the grid at different times (right panel) 
for the De Sitter spacetime. The numerical grid has 28^ points. 
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Figure 2: Convergence of the L2 norm of the Hamiltonian constraint for the De Sitter 
spacetime using two resolutions, one with 14^ and the second with 28^ points on the grid, 
for 20 iterations (left panel) and after 2000 iterations (right panel). 

the Cactus code is given by the time behavior of the Hamiltonian constraint (an output 
of the Cactus code through a thorn called ADMConstraints, namely ham). In the next 
Figure ^ we show the convergence of the L2 norm of the Hamiltonian constraint for 
DeSitter metric, for 20 and 2000 iterations, using two different resolutions on the grid, one 
with 14^ and the second with 28^ points (both grids cover the same region of spacetime, so 
the grid with more points has a smaller value of the finite difference interval Ax). Notice 
that the Hamiltonian constraint for a true solution of the Einstein equations should be 
equal to zero. Finite differencing errors imply that the numerical solution will have 
a non- vanishing value of the Hamiltonian constraint. For a consistent finite difference 
approximation of the Einstein equations, we should expect the Hamiltonian constraint to 
approach zero as the resolution is increased. For a second order approximation, the value 
of the Hamiltonian constraint should go down by a factor of four when the resolution is 
doubled. The figure shows that we have close to second order convergence. 

This was the method to test the convergence of the code in all examples presented in 
this article. We have obtained good second order convergence in all examples as is shown 
in the next figures nr 0, |^, |^, |^ and ^ 
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DeSitter spacetime with cosmological constant 

Convergence of the hamiltonian constraint 

0.0001 



DeSitter spacetime witli cosmological constant 

Convergence of the hamiltonian constraint 

0.0003 




Figure 3: Convergence of the Hamiltonian constraint for the De Sitter spacetime with 
cosmological constant using two resolutions, one with 14^ and the second with 28^ points 
on the grid, for 20 iterations (left panel) and after 2000 iterations (right panel). 

These results were obtained using the BSSN formulation for the numerical evolution, 
but we have obtained similar results using the ADM formulation. As a major conclusion 
we can point out that we found that the BSSN formulation results in stable evolutions. 
For long time evolutions (that is for more iterations than 20, for example 2000 iterations 
in our simulations), the Hamiltonian constraint grows from zero to a constant value, as 
it is obvious by inspecting all our next graphs on the convergence of the Hamiltonian 
constraint (sec below). In our evolutions, one can also observe how the computational 
errors propagate back into the computational grid after being reflected on the boundary 
due to the use of inappropriate boundary conditions (here we used the "flat" boundary 
conditions implemented in the Cactus code). This will be an important problem to solve 
in the future. The De Sitter parameter used in these simulations was a = 1.0 and the 
chosen slicing condition was an "exact" slicing (i.e. use the exact lapse coming from Exact 
thorn). 

Another example of De Sitter cosmology is that of the De Sitter metric with cosmo- 
logical constant, having the form : 

ds^ = -df + e^/^^* [dx^ + dy^ + dz^) . (5) 



where A is the cosmo logical constant. This is an example of an exact solution of the 
Einstein equations (|ID with cosmological constant and without matter. Here we can still 
use the Cactus code, even though it has been written for the case of no cosmological 
constant, by using a simple trick: we have transferred the term with the cosmological 
constant to the right hand side of the Einstein equations. Thus we shall have a non- 
vanishing "matter" term having as components 
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(6) 



(remember that we have G = c = 1). Introducing this in the Cactus code and taking 
A = 1, we have obtained the numerical results presented in the next figures ^ and |[ 
First, we investigated the convergence of the code, using again the time behavior of the 
Hamiltonian constraint. Figure |^ shows this for simulations with 20 and 2000 iterations, 
using again two resolutions, one with 14^ and the second with 28^ points on the grid, with 
< X, 2; < 6 and At = 0.25. We can point out the same conclusions on the convergence 
and the stability of the code as in the case of the De Sitter spacetime. 

Figure ^ presents the evolution in time of the radial metric component grr for 2000 
iterations. It is easy to see here the rapid increase of the radius of the universe (which is 
related to grr\ showing the inflationary nature of a universe modeled by this metric. 



3 Kasner solutions 

We have used the Cactus code for several cosmological solutions of the Einstein equations, 
generically named "Kasner" metrics, included in the class of homogeneous anisotropic 
models, which can emulate some early stages of the Universe. One of these metrics is the 
so-called "Kasner-like" metric (See [17| and |jl9|), which in Cartesian coordinates has the 



form: 
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DeSitter spacetime with cosmological constant 

Evolution of grr for 28^ points on ttie grid 



DeSitter spacetime with csomoiogicai constant 

Values of the gxx at different times on tlie grid 
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time 



Diagonal line of ttie grid 



Figure 4: Behavior of the L2 norm of the radial metric component Qrr (left panel) and 
its evolution in time (right panel) showing the inflationary nature of a universe modeled 
with the De Sitter metric with cosmological constant. 
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Here we have a stress-energy tensor which has all off-diagonal components vanishing: 

/g%# \ 

n .(2 



T- ■ 







q ('-^g)*'' 



V 

This metric forms a one parameter family of solutions of Einstein's equations with a 
perfect stiff fluid. The parameter q is related to the energy density, as is obvious from the 
last equation. The qualitative features of the expansion depend on q in the following way: 
for g > 1/2 the universe expands from a "cigar" singularity; for q = 1/2, the universe 
expands purely transversally from an initial "barrel" singularity; for < g < 1/2 the 
initial singularity is "point-like" and if g < we have a "pancake" singularity. The case 
g = 1/3 corresponds to an isotropic universe with a stiff fluid; the case g = is a region 
of Minkowski spacetime in non-Cartesian coordinates. This family of metrics is "Kasner- 
like" in the sense that the sum of the exponents is equal to one, but the sum of the squares 
is not equal to one except in the cases when g = 0org = 2/3, when we have the vacuum 
case. 

We have investigated the numerical behavior of this metric for different values of the 
parameter g, looking specially at the convergence of the Hamiltonian constraint, starting 
from an initial time t = 1 (since the metric is singular at t = 0), using "fiat" boundary 
conditions and, as in the previous cases, for both 20 and 2000 iterations using the BSSN 
formulation. All convergence tests were done with two resolutions, one with 14^ and the 
second with 28'^ points on the grid. Figure || displays our results for the vacuum case 
(g = 2/3) and 20 iterations. As can be clearly seen the convergence of the Hamiltonian 
constraint is of second order. In the figure we show also the evolution of the L2 norm of 
the Qxx metric component during that time. We show the same plots in figure ^, but this 
time for the value g = — 1 of the Kasner parameter. 

Figure |^ shows the evolution of the radial metric component for the two previous cases, 
pointing out some differences due to the special symmetry of this metric (see above and 
compare with g^x from the previous figures). 

Next we have considered another Kasner type metric, namely the Kasner axisymmetric 
spacetime (jlSl): 

ds^ = ^ + ^ + tdy^ + tdz^ . (9) 

This is an exact solution of the vacuum Einstein equations, explicitly homogeneous, and 
features a cosmological singularity at t = 0. Figure |] shows some of our results for a "long 



(8) 
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Figure 5: Convergence of the Hamiltonian constraint (left panel) and evolution of the L2 
norm of the g^x metric component (right panel) for the Kasner-like spacetime - vacuum 
case {q — 2/3) and 20 iterations. 
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Figure 6: Convergence of the Hamiltonian constraint (left panel) and evolution of the L2 
norm of the Qxx metric component (right panel) for the Kasner-like spacetime for q = —\ 
and 20 iterations. 
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Figure 7: Evolution of the L2 norm of the radial metric component grr for the vacuum 
case {q = 2/3, right panel) and for g = — 1 case (left panel) for the Kasner-like spacetime, 
with 20 iterations in both cases. 
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Figure 8: Convergence of the Hamiltonian constraint (left panel) and the evolution of the 
L2 norm of Qxx (right panel) for the axisymmetric Kasner spacetime using 2000 iterations. 



time" simulation (with 2000 iteration) showing the stability and second order convergence 
of the code. 

A generalization of the Kasner metrics can be introduced with a metric of the form: 



(10) 



where the Kasner parameters pi, p2 and satisfy two relations: 

Pi+P2+P3 = l and p? + P2+P3 = l- (11) 

Restricting ourselves only to two parameters, pi and p2, "we have the following stress- 
energy tensor: 
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where A = pi — p\ + P2 — — P1P2 (note the use of the above first condition on the 
parameters, thus we have = 1 — pi — P2). We have done several simulations for this 
metric. Figure P shows results for 2000 time iterations, using pi = —1/3, p2 = 2/3. 
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Generalised kasner metric 
Convergence test (p1=-1/3, p2=2/3) 



Generalised Kasner metric 

Evolution of the gxx metric component (p1=— 1/3, p2=2/3) 
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Figure 9: Convergence of the Hamiltonian constraint (left panel) and evolution of the L2 
norm of Qxx (right panel) for the generalized Kasner spacetime -pi = —1/3, p2 = 2/3, 
and 2000 iterations. 
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4 Robert son- Walker metric and other examples 



We have studied also other metrics, new or old ones form the Exact thorn. We only 
mention here the Bianchi type I, Godel and Bertotti metrics. A special case is the Kerr 
metric, in two versions, one in Kerr-Schild coordinates and the other one in standard Kerr 
(Cartesian) coordinates. Our results are similar to those discussed above. 

We have further studied the problem of introducing the Robertson- Walker (|16[) metric 
(0) as a generic case for the study of numerical cosmology. Here we use the scale factor 
of the universe R{t) as a variable function through the code, and two initial parameters: 
the initial mass-energy density and the initial scale factor of the universe. Because this 
case has several specific problems and features we shall report on it in a future article in 
preparation. This is intended to be done in connection with the further development of 
the Exact thorn code, where we shall introduce a scalar field coupled with gravitation in 
order to study inflationary models. 
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